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Abstract 

An essential part of describing the damage state and predicting the damage growth m 
a multicracked plate is the accurate calculation of stress intensity factors (SIF’s). Here, a 
methodology and rigorous solution formulation for SIF’s of a multicracked plate, with fully 
interacting cracks, subjected to a far-field arbitrary stress state is presented. The fundamental 
perturbation problem is derived, and the steps needed to formulate the system of singular 
integral equations whose solution gives rise to the evaluation of the SIF’s are identified. This 
analytical derivation and numerical solution are obtained by using intelligent application of 
symbolic computations and automatic FORTRAN generation capabilities (described in the 
second part of this paper). As a result, a symbolic/ FORTRAN package, named SYMFRAC, 
that is capable of providing accurate SIF’s at each crack tip has been developed and validated. 

Nomenclature 

offset angle between inner tips of two parallel cracks 
direction cosines between two local coordinate systems 
strain tensor 

four roots of the characteristic equation 
Poisson’s ratio 

far-held and total stress field, respectively 
components of stress in global coordinate system 
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Vi 

— 

angle defining orientation of local coordinate system 

t,T 

— 

normalized real variables ' 


— 

Fourier transform of Airy stress function with respect to x variable 


— 

half crack length 

a ll 5 ai2i 022 

— 

coefficients of strain-stress relationship (compliant matrix) 

Ui 

— 

auxiliary functions 

*1,*2 

~ 

mode-I and mode-II stress intensity factors 


— 

Fredholm kernels 

Pi 

— 

normal traction at crack surface 


— 

shear traction at crack surface 

s 

— 

Fourier variable 

a 

— 

position vector defining the origin of a local coordinate system 

r jX 1 r iY 

— 

components of the position vector rj 

U, V 

— 

displacement associated with x and y coordinates, respectively 

w 

— 

weight function 

Xj, yj and X, Y 

— 

local and global coordinates 

[A] 

— 

kernel matrix 

C' v c } 

— 

functions of s in Fourier space (i.e., constants in x, y-real space) 

E 

— 

Young’s modulus 

F j( x j, Vi) 

— 

Airy stress function 

Gt)j( t p) 

— 

discrete auxiliary function 

m 

— 

loading vector 


1 Introduction 


In a wide range of materials damage progression is strongly connected with microcrack nucleation 
and growth in solids. Hoagland et al. (1973), Leckie and Hayhurst (1974) and Meyers (1985), 
observed that a large part of applied energy is dissipated because the microcracks extend the 
fracture zone of a larger propagating crack. Similarly, Kachanov (1958) noticed that creep rupture 
is caused by the accelerated growth and localization of microcracks in the material. Therefore, it 
has been concluded that there must be a relationship between continuum damage mechanics and 
fracture mechanics concepts as applied to the problem of a multicracked plate. 

Historically, the solution of a multicracked plane problem is typically solved using one of two 
techniques; the complex potential method ( see Horii and Nemat-Nasser (1985)) or the singular 
integral equation approach (Erdogan (1978)). In this study we have selected the singular integral 
equation technique because of our familiarity with the method and its successful application to a 
number of plane problems; for example the work of Badaliance and Gupta (1976) on isotropic prob- 
lems, Erdogan (1978) and Binienda et al. (1991) on anisotropic and homogeneous problems, and 
Delale (1985) on nonhomogeneous problems involving one or two independent cracks. Although the 
results for multicrack problems obtained using the singular integral equation approach in conjunc- 
tion with the collocation technique are accurate, the required derivation of the system of singular 
integral equations is lengthy and tedious, even for the simplest case. 

One way to avoid this tedious derivation (long hours of searching for possible errors and the 
inevitable “debugging” of a numerical program) is to design and develop a computer code that can 
derive the system of singular integral equations symbolically and can automatically generate the 
associated FORTRAN code required for the solution, by using the collocation technique. With such 
a code a new class of problems can be analyzed more rigorously and efficiently. 
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This subject is covered in two parts; Part I describes the general formulation and solution of 
the stress intensity factors (SIF’s) at the tips of n number of straight microcracks embedded in 
an infinite isotropic plate that is subjected to in-plane loading and Part II (Arnold et al. (1992)), 
describes the symbolic computational aspects and numerical implementation. We expect the sym- 
bolic computational methodology established here and in Part II for n straight cracks in an isotropic 
plate to be easily adapted to an anisotropic material and, in the future, to more complicated crack 
geometries and nonhomogeneous materials. Part I concludes with a comparison to validate the 
accuracy of the results obtained with the present formulation and those reported in the literature, 
whereas other more complex examples, are presented in Part II (Arnold et al. (1992)). 


2 Problem Formulation 

Consider multiple cracks embedded in an infinite isotropic plate (Fig. 1(a)). The plate is under 
a far-field stress denoted by o] x , (in particular a° xx , <t$ y , and where (X,Y) is the global 

coordinate system), and the cracks are defined in their local frames (xj,yj) (Fig. 1(b)). The origin 
of each local frame is defined by the position vector ry, and the orientation of the local frame with 
respect to the global frame is defined by the angle <pj. Each crack is symmetrically situated within 
its own coordinate system and is 2a j long, as shown in Fig. 1(b). 

The general solution formulation can be outlined in four basic steps. The first step is to derive 
the local stress equations for each crack in its respective local coordinate system. This derivation, 
which need be done only once, is achieved by defining the fundamental problem; that is a single crack 
in an infinite isotropic plane (Fig. 1(b)). The fundamental problem is then decomposed into two 
subproblems: the problem of the undamaged plate containing an imaginary crack (Fig. 1(c)) and the 
perturbation problem (Fig. 1(d)) of a plate with a single crack subjected to the appropriate crack- 
surface tractions, which are found from the solution of the complementary undamaged problem. The 
analysis of the perturbation problem leads to singular stresses that govern local crack tip behavior. 

The second step is to formulate the total perturbation stress field for each crack, which includes 
the interaction of all cracks through the summation of the transformed local stresses of all other 
cracks. In the third step of the formulation, the total stress equations are normalized. A set of 
Cauchy type singular integral equations, expressed in terms of unknown auxiliary functions, is 
obtained by subjecting the total perturbation stress equations to the crack-surface traction field at 
each crack location. The fourth and final step of the formulation is to express the SIF s in terms 
of the discrete auxiliary functions, G v j(T p ), evaluated at each crack tip. These discrete auxiliary 
functions are obtained through the implementation of the Lobatto-Chebyshev collocation technique. 
Now let us discuss each step in more detail. 


2.1 Local Stress Formulation 

Consider the fundamental problem (Fig. 1(b)), which is defined as a single crack in an infinite 
isotropic plate, whose solution can be obtained by decomposing it into an undamaged problem 
(Fig. 1(c)) and a perturbation problem (Fig. 1(d)). The essence of this decomposition is that 
the traction forces applied along the crack surface in the perturbation problem are opposite to the 
obtained stress field of the undamaged plate at the particular location of the imaginary crack. As 
a result, this undamaged traction field can be defined in terms of the normal (pj) and shear (qj) 
stress components along the imaginary crack surface: 

Pj(xj) = <r yjV} ( x j,0) ( 1 ) 
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qj(*j) = <r Xi vi(xj, 0) (2) 

« 

where 

a v,v, ( x i » 0) = <T°xx sin Vi + ^y-K cos Vi - <7*y Sin 2<pj (3) 

(xj , 0) = - - — sin 2</?j + a° XY cos 2^ (4) 

The mixed boundary conditions for the perturbation part of the fundamental problem (Fig. 
1(d)) are expressed in terms of stresses 


°y ] v, — Pj{ x j) aQ d — 9j( x j) (®) 

along the crack surface (i.e., y : = 0 and — a } < Xj < aj), and in terms of continuity of displacements 

t> + = v~ and u + = u - (6) 


outside of the crack (i.e., yj = 0 and \xj\ > aj). Here + indicates the value of displacement at a 
point approached from the positive side of the plate, (i.e., y > 0), whereas — indicates the same 
point approached from the negative side of plate, (i.e., y < 0). 

The governing equation for the preceding two-dimensional isotropic plate problem can be ex- 
pressed in terms of the Airy stress function Fj(xj,yj) as 


d 4 F d 4 F d 4 F 

dx 4 dx 2 dy 2 dy 4 


Note: Any stress function that satisfies Eq. (7) and the mixed boundary conditions, Eqs. (5) and 
(6), constitutes the solution to this problem. 

A rigorous solution for this stress function can be obtained by employing the following Fourier 
transformation relation 


y) = ^ J ^ *(*i y) e ,txds 


( 8 ) 


thereby transforming the partial differential Eq. (7) into a fourth-order ordinary differential equation 
with real constant coefficients: 


s 4 $ — 2s 


,cP$ d 4 * 
dy 2 dy 4 


= 0 


(9) 


A known general solution of Eq., (9) is 


$(s,y) = C e* * s (10) 

Differentiating Eq. (10) and substituting the results into Eq. (9), gives the characteristic equation 

p 4 — 2p 2 + 1=0 (11) 

whose roots are pi = p 2 = ~P 3 = — Pa — 1. Therefore, the Airy stress function has the following 
form: 

F(X, y) = jH [(CJ + C' 2 y)e va + (C' + C' 4 y)e-*\ e~ iax ds (12) 

In order to satisfy the physical requirement that the stress function is finite throughout the 
domain of the plate, an alternative form is defined for the upper half plane (for y >0) 

1 f°° . . 

F(x,y+) = — JjC i+ C 2 y)e-^e-*ds (13) 
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( 14 ) 


and for the lower half plane (for y < 0) 

* 

F(x,y-) = ^ f(C 3 + C,y)^t-^is 

Z7T J-oo 

which are automatically bounded at infinity. Note: Constants Cj for j = 1, 2, 3, and 4 are functions 
of the Fourier variable s and are determined by using the local stress continuity conditions at the 
boundaries between the half planes ( y = 0) and by using the perturbation boundary conditions 
subsequent to the determination of the total stresses at each crack location. 

The continuity conditions for local stresses, a yv and <r xy , are identically satisfied, given 


and 


F(x,0 + ) = F(x,(T) 
dF(x, 0+) dF(x,(T) 


dy 


dy 


respectively. Consequently, Eq. (15) gives 

C 3 = C , 

and Eq. (16), upon substitution of Eq. (17), leads to 

Ca = — 2|s|Ci + C 2 


(15) 

(16) 

(17) 

(18) 


Completion of the posed mixed boundary value problem requires the formulation of the local 
stresses and displacements. The local stresses in each half plane are obtained by taking the ap- 
propriate (Timoshenko 1969) second derivatives of the stress function given in Eqs. (13) and (14). 
Thus the stresses associated with the upper half plane are 


+ - J_ 

^ Ir 2 tt 


OO 

J (s 2 C, - 2 |s|C 2 + s 2 yC 2 )e-M» +isx Us 


OO 

J sHCt + c,y)r' M ’+-Us 

— CO 
OO 

< = £ / .(-|*|C, + C 2 - 

— OO 

and those associated with the lower half plane are 

OO 

= I tfC, + 2|s|C, + J J »a)c«*l'-'“l<fa 

Z 7T J 

— 00 

OO 

=-i- J S =(C 3 + C.slell'l*-"'!* 

— CO 
OO 

~ J + c, + \s\yC t )<Mr**is 


(19) 

( 20 ) 
( 21 ) 


( 22 ) 

(23) 

(24) 
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The displacement components (u and v) associated with the x and y directions, respectively, 
are in turn obtained by integration of the appropriate strains, 


«(*.») = J t xx (x,y)di 


(25) 


v(x,y) = J t yy (x,y)dy (26) 

which are related to the local stresses through the utilization of the generalized Hooke’s law: 

Cxr = Qll&xx d" ®120yy (27) 

tyy = a.\lO XX + a\\Oyy (28) 

(29) 

(30) 


where 


and 


a " ~Y - , = T- 


v — 


1 — V 


for plane stress 
for plane strain 


E m = 


E 


for plane stress 
for plane strain 


1 -t/ 2 


(E is Young modulus and u is Poisson’s Ratio). 

As a result the upper half plane displacements are 


oo 

u+ = -ft f h K^i “ 2 \ S \ C * + * a y£a)]e" IWrH- * , <k 

— OO 


and 


+ & f Ts [Ci + <? 2 y] e~ iMy+iax] ds 


OO 

v+ = -%■ J [W<?i + C 2 {\s\y - l)]e-M*+»*l«fc 


OO 

+ f [MIC, + c 2 (\s\y + 1)] e-lWv+^lds 

— OO 

whereas the lower half plane displacements become 

OO 

u" = -%■ / h [(* a Cs + 2M|C 4 + s 2 yC 4 )] e«'lv-.«l ds 

— OO 

OO 0 

+ $* / [C3 + C 4 y]ell J l«-«ld 5 

— OO 

and 

OO 

tr = *£ J [MICa + C 4 (Wy + 1)] e^-^ds 

— OO 

7 [\s\C 3 + C 4 (|s|y - 1 )] cfl-lv—w )ds 


( 31 ) 


(32) 


(33) 


(34) 


It is important to note that the direct application of the stress and displacement boundary 
conditions to the above equations would lead to a set of dual integral equations defined for the 
stresses within the crack domain and for the displacements outside the crack domain. However, in 
this form such a set of dual integral equations is unsolvable for the unknown constants C, and C 2 
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- hence, the introduction of the following auxiliary functions, from which a set of solvable singular 
integral equations can be obtained. The typical form suggested in the literature (Badaliance and 
Gupta 1976; and Erdogan 1978) for these auxiliary functions is, 

/»(*) = - u ‘C 1 ’ 0 )) 

/iW-5l«V,°)-»-(*. 0 )] (36) 

Upon substitution of Eqs. (31-34) and Eqs. (17-18) into Eqs. (35) and (36), respectively, the 
auxiliary functions become 

oo 

/,(*) = 4a„^; / (s 2 Ci - IsICsJe-^ds (37) 


/,(*) = -4<t„ii / |a|a C,e-“'ds (38) 

— OO 

Expressions for the unknown constants C, and C 2 can be determined in terms of the above 
auxiliary functions by inverting the Fourier transformations of Eqs. (37) and (38), where the 
condition that /,(*) and f 3 (t) are nonzero only within the crack region (i.e., -a < t < a) has been 
applied. Otherwise, outside the crack region these functions are zero owing to the displacement 
perturbation boundary conditions; see Eq. (6). Therefore, 


Cl = 



(39) 


Similarly, from Eqs. 
functions: 


a = 


4 a 


li 


-1 

Hi 


f Mt)e ist dt + - r f 2 (t)e iat di 
7-a 5 «/— a 


(40) 


(17) and (18) C 3 and C 4 can also be expressed in terms of the auxiliary 


C 3 


1 

4a n 



(41) 


C 4 = 


-1 

4a» 



(42) 


Substituting expressions for the constants Cj into the local stress equations, for example Eqs. 
(19) to (21), results in the formulation of the set of double integral equations with respect to the 
Fourier variables s (-oo < s < oo) and t {-a < t < a). Integration of this set of double integrals 
with respect to s, will involve expressions having the following recursive form. 


/j = / 0 °° e ps ds where p = p{x,—y,t) 
In = £/n- 1 for n = 2 ’ 3 ’" 


The result of the preceding integration (Gradshteyn and Ryzhik 1980) is 


In 


- — for n = 1,2 3,... 
P n 


(44) 
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For the details of this integration process see Part II (Arnold et al. (1992)). Integrating with respect 
to s will give a set of singular integral equations with respect to t , which are valid for any j th crack 
within its own local coordinate system (xj, yj) : 


T (p) = 


Ua u L, Xi) L» + ((>-*,•)’]» + y] + «! - x,)\ 

r O / . \ o _ a 




y± - (*i - 


(j(i) 

VjVj 


[y] + - x iY ) 2 y] + (<j - x i ) 2 J 

1 / r (4 \(4 \ 2y? 1 

4'7t«h 7-0, V ,2 ^ ,[y> + ( tj - Xj ) 2 ] 2 + y j + ( tj - Xj ) 2 

= 5 ^ - *> ) [[*} + (*?-*,)>]» ‘ yj + «! - ^ 

. , , n »?-(«>- *<)* 1 , 


}*i 


(45) 


(46) 


(47) 


This completes the formulation of the fundamental problem (or local stress state) for the j ih crack. 
Henceforth, the formulation of the multiple crack problem will be addressed. 


2.2 Total Stress Formulation 

The total stress state j{oJ z ) for the j th crack is defined as the local stress state of the j th crack 
( a D pl us the contribution to that stress state of all remaining cracks. This may be represented 
mathematically as 


n — 1 

Vj) ~ Grzfeji Vj) “I" ^rzl x v( x h i/i)? Vv{ x i^ 3/?)1 (48) 

p=l 

for j = l,...,n. Here standard tensor transformation is incorporated, a TZ = fiirflmzGlm and 0i r , fl mz 

are the direction cosines between the (x_,, yj) and (x p , y p ) coordinates with p identifying the remain- 

ing cracks (see Part II). Note: This statement does not imply that the concept of superposition has 
been invoked, since the stress perturbation boundary conditions (see Eqs. (5)) have not yet been 
utilized to determine the unknown auxiliary functions. 

For functional compatibility within Eq. (48), coordinate transformations must be simultaneously 
applied to all remaining (p th ) crack coordinate variables. As a result, the dominant term (i.e., the 
first term of Eq. (48)) possesses a singularity whereas terms within the summation lose their original 
singularities and yet still contribute to the total stress state, as one might expect. The coordinate 
transformation between the (xj, y } ) and (x p , y p ) systems is determined from the following geometric 
relationship (see Fig. 2): 

fjx + Xjcosipj - yjsirupj = r pX + x p cosp p - y p simp p (49) 

rj Y + Xjsimpj + y^cosipj = r pY + x p sirup v + y p cosip p (50) 

where rj X , rj Y are the rectangular components of the j ih crack position vector referred to the global 
coordinate system X — Y, and tpj is the angle of rotation between the global and local systems (see 
Fig. 2). By rearranging Eqs. (49) and (50), one obtains 
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( 51 ) 


x p = for - r pY )sin<p p + fox - r p x)cosif!p 

+xj cos (ip, - (p p ) - yjsm{tpj - <p p ) 

y P = for - r pY )cosp p - fox - r pX )sin<p p 

+Xj sinfoj - <p p ) + yj cos (<p>j - p v ) 


(52) 


2.3 Singular Integral Equations Formulation 

Since the total stress state at each crack location has been formulated, the appropriate stress 
boundary conditions (Eqs. (5)) can be applied. As a result, 2n singular integral equations are 
obtained with 2n unknown auxiliary functions. With the collocation technique, these unknown 
auxiliary functions can be determined in a discrete sense. The collocation technique, however, 
requires that the interval of integration be from -1 to 1 and not from -a to a, as in Eqs. (45-47). In 
order to accomplish this normalization, the following variable substitution, involving the actual j 
crack length a.j, is employed: 

Xj — (lj£ 
tj = djT 

where the £ and r are between -1 and 1. 

The formulation of this system of singular integral equations is complete once the single-value 
conditions for the auxiliary functions, f n j, are chosen. In the case of straight cracks, this single- value 
condition (Erdogan 1978) is: 

J ^fvi( T ) dT = Q ( 55 ) 

where j stands for the j th crack and r/ takes on the value of 1 or 2. 

For example, for the case of an isotropic plate with two cracks, the system of four singular 
integral equations normalized by Eqs. (53 and 54) is obtained by using the cr yy and stress 
components of Eqs. (48) so that 

•< = 7z~i dT + /., keT ”h' dT + /_, ker “ /!i } 

■ = + ker i?/2,<iT + L ker llf22 } 

and 


^ 4a„ | 

n 

fo. 

f h \ dT+ 

1—1 T — £ J 

j ker JJ/ndr + ^ 

f i ker 12 / 12 1 

(58) 

“ 4a„ 1 

n 

I** 

/’ f22 .dr + 

1-1 T — £ J 

f ker jj/ndr + 

^ ker^foj 

(59) 


where ker*? are the Fredholm kernels (Erdogan 1978) and aj,~f,6= 1,2. The above eight kernel 
functions can be represented by four independent relations, given that 6 = — 0 (where 0 = pi — p>\). 
For example, kerj] = ker J}, ker 22 = ker jj, kerj? = ker^, and ker 22 = ker 


(56) 

(57) 


(53) 

(54) 
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For the case when n cracks exist, each singular integral equation would contain more integral 
terms and their associated kernels, thereby giving rise to a system of sihgular equations of dimension 
2n. Note that all of the required kernels for a multiple crack problem will differ from the four kernels 
found for the case of two cracks only by a position vector. Therefore, we can rewrite the total stresses 
(<r^ and <rj v ) for n cracks as follows 


n a xy = 4^17 {/-l kerj/ndr + f± t ker 2 f 12 dr 

+... + /i, ker i/(„_i)idr + fl t ker 2 / (n _ l)2 dr + i /!, ^drj 

n^y = 4^7 {/-I ker s/ll* + f - 1 ker 4 /i idr 

+... + /ij ker 3 /( n _i)idr + /i x ker 4 / (n _ 1)2 dr + £ fl x £$dr} 

(see the appendix for the four independent kernels (ker* for t = 1 ,..., 4 )), and we can associate 
keri = kerj}, ker 2 = ker^, ker 3 = kerjf, and ker 4 = ker^ for the case of two cracks. An 
equivalent association, with differing position vectors, can be done for the case of n cracks. This 
completes the formulation of the system of singular integral equations. 


2.4 Solution for the Stress Intensity Factors 

The integral equations obtained are of the Cauchy type; thus, for sharp cracks the stresses 
and strains will have a square-root singularity and the classic definition of SIF may be used (see 
Badaliance and Gupta 1976; Delale and Erdogan 1977 and 1979; and Delale et al. 1977). Therefore, 
the modes I and II SIF’s for the j th crack are 


*i(l) = lta[2((-l)]* {,•<£((,<>)} 

(62) 

*2<1) = - 1)]> {/t'J.K.O)} 

(63) 

*!(-!) = ( 1 ™,[ _2 ( 1 +f))“{< < ’»(f’°)} 

(64) 

= Hml— 2(1 + 0)1 {,<rj,((,0)} 

(65) 

where the normal and shear stresses, Eqs. (60 and 61), are used. 

It is well known (Erdogan 1978) that the auxiliary functions (/) can be expressed as 
of the unknown bounded functions ( G ) and the known singular weight functions w: 

a product 

/( r ) = G(t)w(t) 

The singular weight function w for a sharp crack is 

(66) 

u»(r) — (t 2 — 1)~ 2 

(67) 

Erdogan (1978) found, for example, that in the case of a Cauchy-type singular integral equation 
(Eqs. (60) and (61)), the dominant part can be expressed in terms of the function G evaluated at 
the tips of the j th crack: 

7 L r- T ( = G + - G ^ T ~ iri ‘ + 0(r) 

(68) 
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where rj is 1 or 2 and O(r) is the higher order term, which in subsequent calculations is neglected. 
The substitution of Eqs. (68) for the dominant parts (the last term in Eqs. (60) and (61)) of the 
normal and shear components of the total stresses in Eqs. (62-65), and subsequent evaluation of 
the limits at the crack tips, results in the redefinition of the SIF’s (normalized with respect to y/ai 
and <Tj p ) expressed in terms of the functions G^j- 


1) = 

4 an v 

(69) 

^(1) = 

4an Glj ( 1 )^° i / ai 

(70) 

fc((-l) = 

7 G2j(-l)\/ a i/ a i 

4a„ 

(71) 

fc^-1) = 

= Gij(—l)\fa.j/a.i 

4 a n 

(72) 


The Lobatto-Chebyshev collocation integration technique is known to provide excellent results 
when dealing with Cauchy-type singular integral equations and so it was used. The unknown 
function G^j is determined at a discrete set of points Tj, T 2 ,...,T m called abscissas. In this way, 
each integral equation is reduced to a set of algebraic equations with unknowns G^j ( T i) 

Gqj ( T m ) , which are the discrete values of the functions G^i] hence its name, a discrete auxiliary 
function. Each of the singular integral equations subjected to the stress boundary conditions (Eqs. 
5), is replaced (see Erdogan 1978) by m - 1 algebraic equations with 2nm unknown parameters: 


1 171 / \ 2 171 

I ■£ t, ” j(Tr p + £ Z + JfemK.) = 4 a, 

V r=l Tt S* 6= 1 r=l 


(73) 


where z,A = a,0,y,S,v = 1,2; p. = fc S; &(* = - 1) are collocation 

points chosen in such a way that ( z / r r ; w r = the weighting constants of the related integration 
formula and Rem = the remainder. For m sufficiently large, Rem can be made as small as necessary 
for the desired accuracy and consequently can be neglected. 

In the Lobatto-Chebyshev method, the abscissas are calculated according to 


(r — 1 )tt , 

Tr - cos — for r = 1, ..., 

m — 1 


m 


with the corresponding weights given by 

7T 


Wi = W m = 


and w T — 


2(m — 1) ' m - 1 

the collocation points are then found by using the formula 


for r = 2, 3, ..., m — 1 


(74) 


(75) 


= cos — — !—■ for z = — 1 (76) 

2tti — 2 

In order to have the complete system of 2nm algebraic equations, the single- value conditions (Eqs. 
(55)) are also expressed by using the collocation technique: 

m 

G v j(T r )w r = 0 (77) 

r=l 

Thus, the resulting system of algebraic equations can be written in the form 
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1-4] {G} = m ■■ (78) 

where [A] is a fully populated 2 nm x 2nm matrix of coefficients and {7£} is the loading function 
vector; see Part II (Arnold et al. (1992)) for details. 

The unknown parameter vector {£?} can be determined through inversion of the [A] matrix; 
thus, 


{G} = [/T'W (79) 

although only the appropriate values (i.e., G^j(±l)) are used to calculate the SIF’s for the j th crack 
(see Eqs. (69-72)). 

Automatic generation of the associated FORTRAN code (see Part II, Arnold et al. (1992)) for 
the evaluation of Eq. (79) completes the development of the solution for any multicrack problem. 
This FORTRAN program was utilized to obtained the following results, which are compared with 
results obtained from other methods existing in the literature. 

3 Numerical Application 

In order to validate the results obtained using the SYMFRAC code, the well-known problem 
of two parallel interacting cracks is considered here. The plate with two cracks of length 2a is 
subjected to a normal stress field (<Tyy) as shown in Fig. 3. The location of the cracks are defined 
by the distance rj between inner crack tips and the offset angle a. 

This problem has been solved for mode I SIF’s and a > 0 by Isida (1976), who used a Laurent 
series expansion that gave less than 2-percent error; by Yokobori et. al. (1971), who used a 
continuously distributed model that gave results within 5-percent error; and by Yijun and Atilla 
(1991), who used an approximate stress field model that gave less than 10-percent error. For the 
case of the collinear cracks (o = 0), an exact solution was obtained by Erdogan (1962) and later by 
Horii and Nemat-Nasser (1985) using the approximate method of pseudo-tractions. 

In Table I the SIF’s at the inner crack tips from Isida (1976), Yokobori et al. (1971), and Yijun 
and Atilla (1991) are compared to those of the present method for the noncollinear case (i.e., a > 0). 
Similarly, in Table II the results for the SIF’s from Horii and Nemat-Nasser (1985) and Erdogan 
(1962) are compared to those of the present singular integral technique for the collinear case (i.e., 
a = 0) at both the inner and outer crack tips. 

Table II shows that the exact results given by Erdogan (1962) and the results from using the 
pseudo-traction method (Horii and Nemat-Nasser 1985) agree exactly with the results given by the 
present singular integral technique. This implies that the error introduced by solving the system of 
singular integral equations with the collocation technique is minimal. 

Also, in the case of noncollinear cracks the SIF’s calculated by the present method are closest to 
the most accurate results of Isida (1976) for the specific cases addressed. We strongly believe that 
the present method is the most accurate because of the minimum number of assumptions required 
for the solution and its favorable comparisons. This fact is further borne out by the observation 
that the most approximate method of Yijun and Atilla (1991) gives the smallest SIF’s for the closest 
crack tip distance, whereas the present method gives the highest values. 
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4 Concluding Remarks 

A rigorous formulation has been presented and validated for calculating the SIF’s of a mult- 
icracked isotropic medium. The size, orientation, and distribution of all cracks were considered 
to be independent parameters of the solution. The formulation was constructed in such a way 
that a symbolic/numeric system (named SYMFRAC, see Part II (Arnold et al. (1992))), capable 
of automatically formulating the required system of singular integral equations for a multicrack 
boundary value problem and generating the associated FORTRAN code, can be developed and uti- 
lized. Extension to anisotropic and/or nonhomogeneous materials, as well as more complex crack 
geometries, would require modifications of the governing equations, boundary /continuity condi- 
tions, constitutive relations, and such. Although the final results would be significantly different, 
the symbolic/numeric procedure developed here would remain the same. Because of symbolic imple 
mentation, the solution of the more complex problems would be a straightforward and automated 
procedure. With this capability, numerous parametric studies could easily be performed to analyze 
the contribution of each parameter on the local stress field as well as the characteristics of the 
damage progression in a material. 

Future work will attempt to generate such solutions and to generate the required data base 
necessary to establish a link between continuum damage mechanics and linear elastic fracture me- 
chanics by homogenizing the effects of interacting distributed cracks into macrodamage models for 
isotropic and anisotropic materials and for more complex crack geometries . These macrodamage 
models can than be efficiently coupled with the finite element method in order to analyze structural 

problems. 
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Appendix - Analytical Form of the Fredholm Kernels 


kerj 


-?{ 


—sin28(p2 + a(stnO)[(a'T — Pi — a(cos0) 2 — (p 2 + a£sin0) 2 ] 
[(p 2 + a£sin0) 2 4- (a'r — pi — a£cos0) 2 ] 2 

2cos29(<i't — pi — a(cos0)(p2 4- a{sind ) 2 
[(P 2 4- a£sin9) 2 4- (a'r — pi — a£cosO ) 2 ] 2 

— sin20(p2 4- a£sinO) — cos29{a'T — p\ — a(cos0 ) 1 
(P 2 + a(sin9 ) 2 -f (a'r — pi — a£co.s0) 2 J 


(Al) 


ker 2 


a' j cos29(p2 + a^sin0)[{atT — pi — a£cos0 ) 2 — (p 2 4- a£smfl) 2 ] 
7r \ [(p 2 4- a£sm0) 2 + (a'r — pi — a£cos0) 2 ] 2 

2sin29(afT — pi — a£co.s0)(p2 4- a£stn#) 2 l 
[(p 2 4- a( i sin0 ) 2 4- (a'r — pi — a£co.s0) 2 ] 2 J 


(A2) 


ker 3 


a' J —2sin20(p2 4- a£sm0) 2 (a'r — pi — a£cos0) 
7r \ [(p 2 4- a£sin9) 2 4- (a'r — pi — a£cos0) 2 ] 2 


co520(p 2 4- a£sm0)[(a/r — pi — a£cos0) 2 — (p 2 4- a£sin0) 2 ] 
[(p 2 4- a£sin0) 2 4- (a'r — pi — a^cos^) 2 ] 2 


2 sin 2 9{p2 4- a£sin0) 4- sin20(a/r — p\ — a£cos9) 1 
(p 2 + a£sin0) 2 4- (a'r — p\ — a£cos8) 2 J 


(A3) 


^ a 1 j —sin29(p2 4- a£.sm#)[(a'r — pi — a£cos#) 2 — (p 2 4- a(iinS) 2 ) 

4 7r \ \{p 2 4- a(sin8) 2 4- (a'r — pi — a£cos6 ) 2 ] 2 

2cos2#(a'r — p\ — a£cos0)(p2 4- a£sin0) 2 
[(p 2 4- a£sin0) 2 4- (a'r — pi — a^cosO) 2 } 2 

(a'r — pi — a£c<>50) 1 

(p 2 4- a(itnS) 2 4- (a'r — pi — a£cos#) 2 J 

where 


Pi = (rj, - r y >)sirup' 4- ( r x - r x .)cosp‘ 
P 2 = ( r y ~ r y /)cos<p' - (r x - r x «)sin<p' 

and 9 = p' — p 


(A4) 


(A5) 

(A6) 
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Table I: Mode-I Normalized SIF’s for Inner Tips of Noncollinear Cracks 


Offset angle, 
a, 

degrees 

r _d 

a 

From literature 

Present 

method 

Yijun and 
Atilla 
(1991) 

Yokobori 
et al. 
(1991) 

Isida* 

(1976) 

45 

1.41 

1.0305 

1.0606 

1.12 

1.1254 

45 

2.82 

1.0135 

1.0533 

1.04 

1.0489 

45 

4.24 

1.0051 

1.0350 

1.02 

1.0254 

38.66 

1.28 

1.0360 

1.0920 

1.13 

1.1317 

38.66 

2.56 

1.0184 

1.0603 

1.05 

1.0551 

38.66 

3.84 

1.0046 

1.0355 

1.03 

1.0302 


* Results taken from graph in literature 


Table II: Mode-I Normalized SIF’s at Inner and Outer Crack Tips of Collinear Cracks 


r. 

From literature 

Present 

method 

0 

a 

Horii and 
Nemat-Nassr 
(1985) 

Erdogan 

(1962) 


Inner 

Outer 

Inner 

Outer 

Inner 

Outer 

0.220 

i 

l 

l 

i 

i 

i 

1.45387 

1.11741 

1.45387 

1.11741 

0.500 

1.2289 

1.0811 

1.22894 

1.08107 

1.22894 

1.08107 

0.857 

1.1333 

1.0579 

1.13326 

1.05786 

1.13326 

1.05786 
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(c) Undamaged plate (d) Perturbation problem 


Fig. 1. Multicracked plate geometry and the method of solution; (a) Multicracked plate, 
(b) Fundamental problem for ] th crack, (c) Undamaged plate, (d) Perturbation 
problem. 
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Fig. 2. Geometric relations between a pair of cracks and their local variables. 
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